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ABSTRACT 



^r I In a series of papers, Nicastro et al. have reported the detection of z > O VII absorption 

features in the spectrum of Mrk 421 obtained with the Chandra Low Energy Transmission Grating 
f— ^ I Spectrometer (LETGS). We evaluate this result in the context of a high quality spectrum of the 

same source obtained with the Reflection Grating Spectrometer (RGS) on XMM-Newton. The 
data comprise over 955 ksec of usable exposure time and more than 2.6 x 10** counts per 50 mA 
at 21. 6A. We concentrate on the spectrally clean region (21.3 < A < 22. 5A) where sharp features 
O I due to the astrophysically abundant O VII may reveal an intervening, warm-hot intergalactic 

4— > ' medium (WHIM). We do not confirm detection of any of the intervening systems claimed to 

jrt ! date. Rather, we detect only three unsurprising, astrophysically expected features down to the 

Log(A^i) ^ 14.6 (3(t) sensitivity level. Each of the two purported WHIM features is rejected 
with a statistical confidence that exceeds that reported for its initial detection. While we can 
^ ' not rule out the existence of fainter, WHIM related features in these spectra, we suggest that 

JH I previous discovery claims were premature. A more recent paper by Williams et al. claims to 

have demonstrated that the RGS data we analyze here do not have the resolution or statistical 
quality required to confirm or deny the LETGS detections. We show that our careful analysis 
resolves the issues encountered by Williams et al. and recovers the full resolution and statistical 
quality of the RGS data. We highlight the differences between our analysis and those published 
by Williams et al. as this may explain our disparate conclusions. 

Subject headings: line: identification — line: profiles — instrumentation: spectrographs — methods: 
data analysis — techniques: spectroscopic — telescopes: XMM-Newton Observatory — Galaxy: halo — 
BL Lacertae objects: individual (Mrk 421) — intergalactic medium — Local Group — diffuse radiation 
— large-scale structure of universe — X-rays: diffuse background — X-rays: ISM — X-rays: individual 
(Mrk 421) 

1. Introduction 
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servatories (Canizares et al. 2000; Brinkman et al. 
2000; den Herder et al. 2001), it has become fea- 
sible to undertake an exploratory search for the 
large amount of baryonic matter that may be 
contained in a highly ionized phase of the lo- 
cal Intergalactic Medium. The agreement of the 
baryon density at high redshift as measured in 
the Lya forest, compared to the density predicted 
from Big Bang Nucleosynthesis and the light el- 
ement abundances (Cowie et al. 1995; Buries & 
Tytler 1997, 1998) strongly suggests that such a 
medium should exist. Independently, large scale 
coupled dark matter/hydrodynamics simulations 
have shown that an early, largely neutral IGM will 
progressively become more highly ionized, until at 
the present day it is essentially undetectable at op- 
tical/UV wavelengths (Cen & Ostriker 1999; Croft 
et al. 2001). The calculations indicate that a ma- 
jor fraction of the baryons at small redshift could 
reside in this diffuse, warm, unvirialized phase of 
the IGM, and this is consistent with the fact that 
the local baryon density inferred from a census of 
stars and gas in virialized structures falls short of 
the predicted value by up to 50%. 

Currently, the most promising technique for de- 
tecting and characterizing the medium is high res- 
olution soft X-ray absorption spectroscopy of the 
low-Z elements towards bright extragalactic con- 
tinuum sources, which has now been attempted 
with a variety of instruments towards a number 
of suitably bright, spectrally featureless objects 
(Fang & Canizares 2000; Fang et al. 2001, 2002a,b; 
Nicastro et al. 2002; Rasmussen et al. 2003a, b; 
Mathur et al. 2003; Cagnoni et al. 2004; Ravasio 
et al. 2005; Nicastro et al. 2005b; Barcons et al. 
2005). The search has naturally focused on the 
K shell resonance lines of H- and He-like oxygen 
(O VII ls2 i^o - ls2p iPi (w), 21.602 A, and 
O VIII Is — 2p, 18.969 A), since oxygen has high 
abundance, and the O K band is relatively clean. 
Thus far, these observations have not produced 
an unambiguous detection of highly ionized met- 
als except at zero redshift, where resonance ab- 
sorption in C, O, and Ne is associated with gas 
in and around the Milky Way Galaxy, and pos- 
sibly in a tenuous intragroup medium in the Lo- 
cal group. This state of affairs is not surprising, 
given the predicted distribution of column densi- 
ties (Fang et al. 2002a; Chen et al. 2003), the fact 
that only a relatively small redshift path has been 



surveyed, and that the wavelength resolution of 
the grating spectrometers on Chandra andXMM- 
Newton is roughly an order of magnitude too poor 
to provide adequate sensitivity for the expected 
equivalent widths. 

The most ambitious and most promising search 
has been conducted using a very deep spectrum of 
the blazar Mrk 421, obtained with the Low Energy 
Transmission Grating Spectrometer [LETGS) on 
Chandra (using both the ACIS-S and HRC-S cam- 
eras) , collecting data taken during times when the 
source was undergoing an outburst (Nicastro et al. 
2005a, b). The most surprising feature in this spec- 
trum is the presence of what appear to be faint 
O VII w (21.602 A) resonance absorption lines, at 
redshifts z « 0.011 and z w 0.027 (the redshift of 
Mrk 421 (Ulrich et al. 1975) is z w 0.0308), and ad- 
ditional absorption, though much weaker, in other 
transitions at approximately the same redshifts. 
These reported detections have attracted substan- 
tial attention, since, if correct, they would repre- 
sent the first discovery of the hottest phase of the 
WHIM, long predicted by cosmological N-body 
simulations. 

We have analyzed an even deeper spectrum 
based on data obtained with the Reflection Grat- 
ing Spectrometer {RGS) on XMM-Newton, com- 
prising nearly a million seconds of exposure time, 
in a spectrometer with approximately 3 times the 
effective area of LETG/HRC-S (5-8 times the 
effective area of LETG/AGIS-S) in the relevant 
spectral band, at comparable wavelength resolu- 
tion. We do not confirm the detectiion of the 
z > O VII absorption features reported by Nicas- 
tro et al. (2005a,b). 

In a recent paper, Williams et al. (2006) re- 
port on the analysis of a portion of these same 
RGS observations and claim to have shown that 
the data are not of sufficient quality to confirm 
or deny the Chandra detections. On the contrary: 
We show that a careful reduction of the RGS data 
does not suffer from the serious limitations cited 
by them, and that the RGS instrument may have 
been incorrectly blamed for the shortcomings of 
their attempted analysis. 

In the following, we describe the RGS dataset 
on Mrk 421 and we give a detailed description of 
our data analysis procedures. We then describe 
the search for faint, redshifted discrete absorption 
lines in the RGS spectrum, and show that no sig- 



nificant redshifted absorption is seen in the O K 
band. We quantify our result that absorption hncs 
at the contrast seen in LETGS can be ruled out. 
We compare our analysis with that of Williams 
et al. (2006) and illustrate what we believe the dif- 
ferences are. We conclude with a summary, and 
briefly put our finding in the context of the search 
for the highly ionized IGM. 

2. Observations 

Already in its sixth year in orbit, XMM-Newton 
has pointed toward Mrk 421 many times for vari- 
ous purposes. Our strategy to maximize our sen- 
sitivity for detecting faint spectral features logi- 
cally requires combining as much data as possi- 
ble, drawing from the multiple observations that 
are available. Table 1 summarizes the 33 XMM- 
Newton pointings (out of 36 currently available) 
that we used in our analysis. We combined data 
from these observation data files (ODF) into 14 
different data sets. After rejecting periods of 
high background, the remaining "good" integra- 
tion time makes up roughly 93% of the total (ex- 
ceeding 1 Ms), at 955 ks. Three ODFs that were 
neglected would have contributed only 10 ks in 
additional exposure. 

The data sets that were analyzed are summa- 
rized in Table 2. Essential data products (spec- 
tra and response matrices) for each data set were 
prepared using filtering parameters specific to the 
content of the constituent ODFs. A description of 
the filtering specifications is included below. 

The usual data analysis paradigm that applies 
to our approach is one where a data set's integra- 
tion time occurs over a period where the source 
and all instrument characteristics are assumed to 
be static. The reality is substantially different, 
particularly for long observations. Even in a sin- 
gle observation, many parameters of an observa- 
tion change, including the source's spectrum, the 
satellite's attitude, and the number and location of 
problematic detector areas that are identified and 
subsequently rejected from analysis. Variation of 
any of these quantities can lead to subtle compli- 
cations and may frustrate simple interpretation of 
weak features in the data. The proper treatment 
of these effects is detailed below. 



3. Data Analysis 

The RGS branch of the XMM-Newton Science 
Analysis System (SAS^) provides software filtering 
for X-ray events detected in the readout CCDs. It 
is designed to provide nominal data products es- 
sential for spectral modeling (spectra, background 
samples and response matrices) under the assump- 
tion that the X~ray event data can be converted 
into a form where the instrumental signature of 
each event's origin has been removed. While 
the response matrix generator can provide proper 
compensation for physical QE variation across 
the CCD array, it specifically does not compen- 
sate for QE variations due to local charge trans- 
port anomalies across the 18 detector readouts per 
RGS. As the CCDs incur more radiation damage, 
this becomes a limitation of the existing analysis 
approach. Thus, while overall charge transport 
characteristics are corrected for, the event "pulse- 
invariant" conversion model does not include de- 
tail on the pixel scale. Currently, only four charge 
transfer inefficiency (CTI) parameters per readout 
node are used to perform this correction. 

We find that the majority of faint, systematics- 
induccd features arise from three mechanisms: 

1. Localized gain anomalies (attributed to ra- 
diation damage induced CTI variations). 
These are identified and corrected for in 
hardware coordinates. In general these are 
limited to a small fraction (<1%) of CCD 
area, particularly after camera cooldown was 
performed circa rev. 0530. 

2. Transient, high duty-cycle pixel reads. 
These can evade identification in long ex- 
posure times, where the hot pixel finder will 
succeed in finding persistent pixels even at 
a low duty cycle. 

3. Crosstalk pixels - pickup of synchronously 
sampled analog signal of high dark current 
pixels. The CCDs are each read out of two 
output amplifiers and "mirror images" of 
cosmic rays are seen in the electronic image 
of the other readout, at a ^1-2% crosstalk 
coupling. 

4. Changes in source spectrum in the presence 
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Table 1 

MRK 421 OBSERVATIONS FOR ABSORPTION ANALYSIS. 



XMM-Newton 


Obs. Start 


Obs. On 


Target 


Obs. ID 


Date (UT) 


Time [s] 


Dataset 


0099280101 


2000-05-25 


66497 


0084 


0099280201 


2000-11-01 


40115 


0165 


0099280301 


2000-11-13 


49811 


0171 


0099280401 


2000-11-14 


43010 


0171 


0099280501 


2000-11-13 


21206 


0171 


0099280601 


2000-11-15 


20213 


0171 


0136540101 


2001-05-08 


39007 


0259 


0136540201 


2001-05-08 


9816 


0259 


0153950601 


2002-05-04 


39727 


0440 


0153950701 


2002-05-05 


19982 


0440 


0153950801 


2002-05-05 


21671 


0440 


0136540301 


2002-11-04 


23913 


0532 


013654040P 


2002-11-04 


23917 


0532 


013654050P 


2002-11-04 


22914 


0532 


013654060P 


2002-11-04 


22917 


0532 


0155555501^ 


2002-11-05 


37765 


0532 


0136540701'" 


2002-11-14 


71520 


0537 


0136540801'" 


2002-11-14 


11415 


0537 


0136540901*" 


2002-11-15 


11420 


0537 


0136541001 


2002-12-01 


71118 


0546 


0136541101 


2002-12-02 


11413 


0546 


0136541201 


2002-12-02 


11415 


0546 


0158970101 


2003-06-01 


47538 


0637 


0158970201 


2003-06-02 


8963 


0637 


0158970701 


2003-06-07 


8055 


0640 


0158970801 


2003-06-07 


12805 


0640 


0158970901 


2003-06-08 


10752 


0640 


0158971001 


2003-06-08 


12800 


0640 


0150498701 


2003-11-14 


48917 


0720 


0162960101 


2003-12-10 


39889 


0733 


0158971201 


2004-05-06 


66141 


0807 


0153951201 


2005-11-07 


10017 


108x 


0158971301 


2005-11-09 


60015 


108x 


•^RFC (RGS2) 


Cooldown 






''RFC (RGSl) 


Cooldown 







Table 2 
Spectral data sets produced. 



Dataset 


GTI 


N a 
AA 


Aflgob 

R 


A^Z-go"^ 


S^ 


xl^ 


0084 


63.52 


1670 


0.29 


0.004 


0.001 


0.912 


0165 


36.34 


426 


0.33 


0.010 


0.003 


0.900 


0171 


133.03 


4383 


0.25 


0.241 


0.061 


1.197 


0259 


40.84 


1007 


0.23 


0.010 


0.002 


1.042 


0440 


80.86 


862 


0.62 


0.107 


0.067 


0.979 


0532 


92.31 


3044 


0.52 


0.050 


0.026 


1.331 


0537 


91.04 


3028 


0.48 


0.169 


0.081 


1.028 


0546 


93.66 


1722 


0.29 


0.168 


0.049 


1.000 


0637 


62.95 


1276 


0.26 


0.104 


0.027 


1.051 


0640 


49.02 


584 


0.18 


0.517 


0.096 


1.194 


0720 


48.80 


2094 


0.28 


0.024 


0.007 


1.054 


0733 


27.78 


572 


0.17 


0.012 


0.002 


0.923 


0807 


65.98 


2632 


0.28 


0.006 


0.002 


1.222 


108x 


69.82 


2888 


0.23 


0.102 


0.023 


0.770 


Total 


955.90 


26188 








1.040 



'^counts per 50 mA spectral range at 21. 6A 

''90% widths of the relative RGS count rate distribution, 

ofF-axis angle distribution, and their product S = Aipgo x 
A-R90 

R 

=fits to the AA=1.2A band, with 7.4 mA bins 



of finite spacecraft drift. Systematics are in- 
troduced in regions of non-continuous detec- 
tor coverage. The effect is present only if 
both spectral variation and spacecraft drifts 
occur within an observation. It is analogous 
to the effect of adding spectra from multi- 
ple offset pointings where the source varied 
between pointings (c/. § 5.3). 

The first three mechanisms listed above have been 
observed with varying degrees of significance, and 
each provides a mediating process for redistribut- 
ing nominal X-ray induced CCD events out of the 
event extraction pulseheight window (thereby in- 
ducing faint systematics that can resemble absorp- 
tion features). The fourth mechanism can affect 
the interpretation of the data in the context of the 
common approach to analysis^. Three columns 
in Table 2 address some of these effects ( ^^" , 
A09O & S). These correspond to the 90% distri- 
bution widths in the relative countratc and the 
pointing variation (in arcminutes) and their prod- 
uct {S = A0go X ''^"' ), respectively. Contribution 
to systematics by this mechanism may scale with 
a data set's susceptibility S and the number of 
affected regions should scale with the number of 
detector regions excluded from the analysis. 

Our production of data sets that are minimally 
impacted by systematics focuses on identification 
and rejection of problematic detector regions and 
is guided by the following principles: 

1. The effective exposure time of each data 
set should be sufficient to generate adequate 
statistics for a quiescent thresholded hot 
pixel map and median offset map for each 
CCD. Conversely, the instrumental detail 
for each data set should be specific to its 
epoch and not averaged over other obser- 
vations where detector characteristics may 
have changed. We have chosen to com- 
bine ODF data that span approximately full 
XMM-Newton revolutions (^48 hours). 

2. An iterative approach is adopted for rejec- 



^ Consider co-adding raw spectral data from two phases 
{e.g., dim and bright) of a source spectrum, where gaps 
in spectral coverage have also moved between the phases. 
When a single exposure map is used in the analysis, bi- 
modal features are induced in the residuals to any spectral 
fit, whose locations correspond to the gap's position in each 
of the two phases. 



tion of problematic detector areas. Auto- 
matic detection and rejection of some re- 
gions is performed by simultaneously in- 
specting both a median offset CCD frame 
and a thresholded hot pixel map for each 
CCD analyzed. A second screening is per- 
formed by manual inspection of the sur- 
viving X-ray events, arranged by hard- 
ware and pulse-invariant event parameters. 
Non-statistical features are readily identi- 
fied when superior counting statistics are 
available. An example for the iterative ap- 
proach to data screening is provided in the 
Appendix. 

While we used custom software to prepare the 
datasets for spectral analysis, we emphasize that 
if used thoughtfully, any software can be used to 
achieve these results. In the Appendix we pro- 
vide a demonstration that the XMM-Newton SAS 
can be used (interactively and iteratively) to pro- 
duce an equivalent spectrum. We cannot guaran- 
tee this if one were to start with the data products 
automatically generated by the SAS pipeline pro- 
cessing system (PPS). This limitation is consistent 
with the purpose of the PPS products: They are 
intended to provide the user with a "quick-look" 
assessment of the ODF's contents, and no attempt 
is made in their preparation to reduce systematics. 

Each data set generated was then analyzed in 
parallel within XSPEC (Arnaud 1996) as a sep- 
arate spectrum paired with its specific response 
matrix and backlight emission model. A total of 
10 isolated channels (out of a total 2268) stood 
out by contributing A^^ > 10 and were excluded 
from the fit. Indeed, 8 of these 10 channels reside 
in the 6 datasets with the susceptibility parameter 
S > 0.025' (c/. Tab. 2), consistent with some of 
our initial considerations on systematic errors. 

The spectra of all the data sets included in the 
joint fit are combined only for plotting purposes 
via XSPEC's setplot group command. The re- 
sulting absorption spectrum is given in Figure 1, 
where each spectrum was divided through by its 
specific folded powerlaw continuum model prior to 
averaging across data sets. In the following section 
we discuss the quantitative analysis of this com- 
posite spectrum representing 955 ks of integration 
time toward Mrk 421. 



4. Analysis of the 21.2 22. 5A spectrum 
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Fig. 1. — The deep (955 ksec) absorption spec- 
trum toward Mrk 421 (top panel). Locations 
for intervening absorption features previously re- 
ported by Nicastro et al. (2005b) are indicated 
with parenthesized redshifts. The Mrk 421 spec- 
trum is given alongside a line-rich coronal spec- 
trum (bottom panel) of comparable integration 
time and observation count for comparison. Both 
plots use the instrument's blind wavelength scale 
for their constituent datasets. The O VII w emis- 
sion line profile in the spectrum of AB Dor pro- 
vides a model- and calibration-independent esti- 
mate for the effective wavelength response in spec- 
tra composed of many different datasets. 



4.1. Continuum modeling; the discrete ab- 
sorption at redshift zero 

The absorption spectrum in the range where 
the O VII w line should appear for the redshift 
range z = — 0.030 was generated as described 
above. No additional splines to the instrument 
response were applied, so any uncorrected, differ- 
ential calibration errors over this spectral range 
will appear in these residuals. We neglected any 
absorption by neutral gas in our Galaxy in the 
fitting, since our analysis spans a narrow wave- 
length band. The foreground column density is 
small (Ah = 1.61 x 10^" cm~^, Lockman & Savage 
1995), and the transmission T of the ISM varies 
by 2% across the chosen range, but can be approx- 
imated locally by a power law T ~ ^-0-38 -^^jthin 
0.02 % accuracy. The most prominent absorption 
features (A = 21.6, 22.0 A) were fitted at the same 
time, but the wavelength and strength were fixed 
across all data sets because the absorber is consid- 
ered to be independent of the backlight spectrum. 
Our method therefore reduces any differences in 
the continuum between data sets. 

Next, absorption parameters for each line fea- 
ture were further determined by fitting only in a 
narrow spectral range centered on each feature. 
The powerlaw indices of each data set (determined 
by fitting over 21.3 < A < 22.5 A) were held fixed 
and only their normalizations were allowed to vary. 
The fitting ranges were (21.4 < A < 21.8 A), 
(21.9 < A < 22.1 A) and (22.2 < A < 22.4 A), 
respectively, for the three detected features listed 
in Table 3. 

To fit absorption line profiles and to estimate 
equivalent widths, we generated response matri- 
ces with energy bin densities that oversample the 
RGS spectral resolution by a factor of about 5 
(AA ^ A/1800), sufficient to model narrow fea- 
tures in the spectra. For the intrinsic line profiles 
we used the Voigt profile, evaluated appropriately 
for a discrete model wavelength grid. 

The fit parameters providing the best fits (for 
given temperature, atomic mass and Doppler 
parameter) are given in the oscillator strength- 
column density product (/A^) and the equivalent 
width Wx for the line (again, computed on the 
spectral model grid). 



Table 3 
Joint fitting results and line detection limits for datasets in Table 2 



Feature X„ 



AX"" 



ID^ 



Wa^ 



S/N'^ 



LogNid X^(d.o.f.) 



21.595 


+ 3.2 
- 3.0 


VII 


= 


13.6i- 


25.6(20.4) 


1 ^ 7fi'?+0031 
I0./D^)_o.o46 


0.951(726) 


22.022 


+ 10.0 
-10.0 


VI 


0.0004 ± .0005 


3.3t- 


6.3 (4.2) 


15.22 t'of. 


0.902(332) 


22.290 


+21.1 
-19.9 


(0 VII) 


0.0318 ±.0010 


2.it}:^ 


4.0 (2.6) 


14.93 t^'il 


1.190(347) 









Non-detect 


;ions 








AA'' 


ID'^ 


z° 


Wa" 


[Log Nil f 


xUd.o.L) 


21.85 
22.20 
22.32 


+20 

+20 
-20 
+20 
-20 


(0 VII) 
(0 VII) 
(0 VII) 


0-028t-Z 


+^-^ 

^ -1.4 

+2-0 
"-12 

^ -0.5 


14.78(15.01) 
14.81(15.06) 
14.93(14.93) 


1.025(346) 
1.229(347) 
1.182(348) 



''Quoted uncertainties are for 90% confidence limits (Ax^ = 2.706). Wavelength uncertainties and 
equivalent widths are given in mA. 

'^Parenthesized entries are tentative. O VII and O VI represent the Ka blends A21.602A and A22.02A, 
respectively. 

"^Significance of each detection. The first number given is yAx^ in the well-sampled continuum limit. 
Parenthesized values are for narrow spectral range fitting: AA =:0.4A for Feature 1, AA =0.2A for other 
features. 

'^Column densities corresponding to each feature. Feature 1 was fit best when using a turbulent velocity 
parameter v — 330kms~^; the other two lines were fit using a fixed v = lOOkms"^. The effective oscillator 
strengths assumed for O VII and O VI were 0.695 and 0.525, respectively. 

''Quantities were taken from the detections of Nicastro et al. (2005b). 

'90% confidence upper limits to absorption column density. The first number corresponds to fitting 
the line feature at the input wavelength tablulated by Nicastro et al. (2005b); the parenthesized value 
corresponds to allowing the absorber wavelength to vary freely within the 40 mA band centered on Xinput , 
the wavelength uncertainty they assumed for Chandra LETGS. 



Three apparent absorption features stand out 
visually: at wavelengths A w 21.60, 22.02, 22.29 A. 
The results of formally fitting these features are 
summarized in the first three entries of Table 3. 

The strongest line is identified with the the 
O VII w transition at z = (21.602 A). The best 
fit wavelength is off by 7 mA (or 100 kms^^) from 
the expected value and formally outside of the 
measured 90 % confidence limit range, by 4 mA; 
however this mismatch is smaller than the overall 
wavelength scale zero point uncertainty. Therefore 
we assume that it is the 21.602 A line feature at 
z = 0, and it provides a local wavelength fiducial 
for the remaining features in the spectrum. 

The next strongest line is coincident with the 
dominant inner shell excitation doublet of O VI, 
also at 2; = 0, which has a laboratory wavelength of 
22.0194±0.0016 A (Schmidt et al. 2004). This line 
is also quite prominent in Figure 1. The measured 
strength of this line {W\ ~ 3.3 mA) is marginally 
consistent with the Chandra LETGS measniement 
of 2.4 ± 0.9 mA (WilHams et al. 2005) but the 
inferred column density raises some concern, be- 
cause it is large compared to the column density 
in O VI resolved by FUSE (Savage et al. 2005). 
The narrow feature immediately to the right of the 
O VI inner shell feature appears to be narrower 
than the instrumental resolution. The source of 
this fluctuation could be instrumental, but an ex- 
haustive search for its origin has not been per- 
formed. Unfortunately the wavelength range is 
only covered by RGSl^ and an independent cross 
check on this feature is not available. 

Finally, the weakest, ^ 2ct feature close to 
22.29 A is suggestive of absorption by O VII 
21.602 A in the host galaxy (0.031 < z < 0.033). 
The 90% confidence range for the feature redshift 
is marginally consistent with a currently accepted 
redshift of z = 0.0308 for Mrk 421 (Ulrich et al. 
1975). They did not give an accuracy estimate for 
their redshift but an inspection of their published 
plot suggests it may be of the order of 0.001. 

4.2. Weak line feature search 

To identify any other possible spectral features 
in the wavelength range of interest, we have per- 



•^See Appendix B for a study of typical systcmatics intrin- 
sic to the data, performed by comparing differences in the 
recorded spectrum between the two RGS models. 



formed a search by fitting for lines within localized 
spectral ranges. The fitting ranges were 4 resolu- 
tion elements in width (0.2 A) and were centered 
on the hypothetical feature wavelength. Allowing 
only the normalizations of the local continua to 
vary and with the powcrlaw indices fixed, the nar- 
row spectral range was fitted. Then a line feature 
(either in absorption or emission) was introduced 
to the model and its best fit amplitude was deter- 
mined, along with the change in x^- This process 
was performed on a fine, uniformly spaced grid 
of wavelengths within the search region. Results 
of this process are given in Fig. 2. The search 
is effectively a localized test of the null hypoth- 
esis (no line), performed over the spectral range. 
While some features appear to be fit with signifi- 
cances of 2cr, the overall distribution of these "de- 
tections" is well behaved (Fig. 3) and the number 
of 2(T detections is not greater than expected, ac- 
cording to the x^ and F distributions. The typical 
null hypothesis probability for obtaining our value 
for F is 8.5%. We conclude that at the working 
sensitivity level of the blind search, the data are 
consistent with no excess in features above the 2a 
level (1.9 mA; Nqvii ^7 x lOi'^cm'^) that may 
trace intervening WHIM toward Mrk 421. 

5. Comparison to Published Results 

5.1. Comparison to Williams et al. (2005) 

We immediately notice that our measured W\ 
values are significantly larger (>30%) than those 
tabulated based on the Chandra results. We in- 
vestigated reasons for this, and found that the 
discrepancy is most likely due to the combined ef- 
fect of several factors. First, it should be noted 
that W\ values derived from order-unsorted data 
(LETG/HRC-S) can depend strongly on proper- 
ties of the underlying continuum model and the 
assumed column density. When we fit the LETGS 
data ourselves (Kaastra et al. 2006), we obtain W\ 
values for the order sorted (LETG/ACIS-S) obser- 
vations that exceed the published (Williams et al. 
2005) values by roughly 20%. Moreover, when we 
attempt to estimate Wx values graphically from 
their Figure 1(c) we are able to reproduce their 
tabulated values precisely. We found this surpris- 
ing because the presence of spectral contamination 
in their composite data set (estimated at about 
20%) should reduce the apparent Wx value cor- 




feature input wavelength [A] 

Fig. 2. — Results of the search for narrow, un- 
resolved features in the O VII region of the spec- 
trum. The positive detections identified in Table 3 
have been included in the underlying model. The 
top plot shows yAx^ (in units of a) between the 
best fit local continuum level and adding in a line- 
like feature at the center of the data band where 
the fitting was performed. The bottom plot gives 
the best fit parameter fNi which encodes the po- 
larity and equivalent width of the best fit feature. 
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Fig. 3. — The distribution of "detections" in the 
narrow feature search results (c/. Fig. 2). The 
distribution appears to be consistent with what 
is expected, based on the y^ and F distributions 
(the overplotted curve). There is no indication for 
additional intervening absorption systems at this 
effective sensitivity level {a corresponds to 5W\ ~ 
0.95 mA). 



respondingly. We therefore believe that contam- 
inating contributions to the continuum level had 
not been properly subtracted prior to their W\ es- 
timation, thereby resulting in systematically low 
equivalent widths. Because the discussion of their 
data analysis refers to Nicastro et al. (2005b) for a 
thorough description, we infer that the systematic 
underestimation affects both of these papers. 

The presence of finite spectral contamina- 
tion by scattering further decreases apparent 
equivalent widths. While integrated contri- 
butions of this sort are considered small, we 
note that off-diagonal terms are neglected in 
off-the-shelf Chandra matrices: Response val- 
ues (e.^., http://asc.harvard.edu/cal/Links- 
/Letg/User/HrcjQE/ea_index.html) are identi- 
cally zero for off-diagonal elements where |mAA| > 
75 mA. 

On the i?G5' side, it is certainly conceivable that 
estimated W\ values may be systematically high 
or low by a relatively small (~10%) amount. Pre- 
fiight calibration activities included illuminating 
individual RGS gratings to measure their scatter 
distribution over various angular scales. Interpre- 
tation of the raw data naturally placed upper lim- 
its to true scatter contributions, but improper ac- 
counting for of any spectral contamination in the 
source or scatter in the beamline each affect the 
apparent scatter amplitude off of the gratings. A 
modest overestimation in the grating scatter am- 
plitude in the RGS physical model also lead to 
overestimates in W\ measures, according to our 
method, by approximately the same amount. 

It is evident, therefore, that the z—Q O VII fea- 
ture in the spectrum is reported inconsistently in 
units of mA. Unless otherwise noted in the follow- 
ing discussion, we disregard this fact and use tab- 
ulated entries for absorption line strengths {W\) 
taken at face value. 

5.2. Comparison to Nicastro et al. (2005b) 

As we have already shown above, a blind search 
of the RGS spectrum for narrow absorption lines 
in the 21.3 — 22.5 A wavelength range yields no ev- 
idence for the presence of lines that could indicate 
absorption by intervening WHIM. Furthermore, 
directly fitting for the absorption features identi- 
fied in the Chandra LETCS spectrum reported by 
Nicastro et al. (2005b) yielded only null results (c/. 
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the three non-detections tabulated in Table 3). 

We are left to compare the XMM-Newton RGS 
absorption spectrum toward Mrk 421 wc obtained 
to the folded spectrum that is consistent with the 
reported intervening absorbers. This is given in 
Figure 4. To generate this comparison we took 
the feature equivalent widths tabulated in the two 
LETGS papers (Nicastro ct al. 2005b; Wilhams 
et al. 2005) together with a quantitative estimate 
for the z = 0.033 feature which was not tabulated 
in either paper, to produce a pattern of absorp- 
tion features. The pattern amplitude was then 
scaled to obtain agreement between the data and 
the folded model for the O VII w (z = 0) feature 
as given in Table 3. 

The comparison between the RGS data and 
the synthesized LETGS model (folded through the 
RGS response) is striking. Overall, two features 
are in reasonable agreement, while three features 
are not. All three features that are inconsistent are 
the purported intervening O VII w absorption sys- 
tems at redshifts z «0.011, 0.027 & 0.033. Their 
presence, at the contrast seen in the LETGS spec- 
trum is ruled out at the 3.5 ct, 2.8 ct and 5.2 ct lev- 
els, respectively'*. We can conclude that the fea- 
tures in the Chandra data as reported by Nicastro 
et al. (2005b) are in disagreement with the results 
of our analysis of the RGS data. 

5.3. Comparison to Williams et al. (2006) 

Williams et al. (2006) searched for the z> 
absorption lines in a subset of RGS data. They 
attributed their inability to confirm the interven- 
ing WHIM features discovered by Nicastro et al. 
(2005b) to a host of instrumental problems intrin- 
sic to the RGS, and declared their non-detection 
fully consistent with the Chandra measurement. 
We have demonstrated that a careful analysis is 
not in agreement with this conclusion. We have 
successfully used the XMM-Newton RGS to pro- 
duce and analyze data sets that are minimally af- 
fected by instrumental systematics. Our strong 
exclusion of only the redshifted features detected 
by Chandra makes the claim for the discovery of 
the WHIM dubious. 

By studying the analysis approach described in 
Williams et al. (2006), we have ascertained that 



they introduced certain problems in the method 
for co-adding the 14 observations^. For a source 
as variable as Mrk 421, where the spectra are ac- 
cumulated from a large number of pointings char- 
acterized by small offsets, certain systematics will 
naturally be introduced into the data wherever 
bad detector regions exist®. The key lies in un- 
derstanding how non-continuous detector cover- 
age can be used to still yield results that are not 
riddled with artifacts from that coverage. Because 
the countrates in the RGS varied by as much as a 
factor of 4 for the observations they used, a sin- 
gle bad column that was properly ignored in the 
data pipeline (but improperly compensated for in 
modeling) can conceivably introduce multi-modal 
features with equivalent widths as large as 6 niA''. 
This is a factor of 10 or so greater than the 1 ct 
counting statistics limit. 

While we accumulated counting statistics to 
a different level (26000 vs. their 12500 cts per 
50 mA), this fact alone should have had an in- 
consequential impact on the overall sensitivity to 
spectral features of the significance detected by 
Chandra LETGS. We suggest that Williams et al. 
(2006) did not exercise proper care in collecting up 
large quantities of data for the purpose of measur- 
ing faint spectral features. 

Irrespective of our mutually opposing conclu- 
sions, secondary reasons for obtaining different 
quality spectra using the RGS exist, and we list 
them here: 

1. Our data were not added together, but fit- 
ting residuals were averaged for display pur- 
poses. Each dataset was analyzed with its 
corresponding response matrix, as described 
above. This naturally led to nearly complete 
spectral coverage over the range of interest. 

2. Large countrate fluctuations ('--^5%) from 
channel to channel are seen in their spec- 
trum (their Fig. 1), and are suggestive of 
aliasing problems. Similar fluctuations are 



*These exclusions are based on Wx face value tabulations 
and not on the line strengths reflected in Fig. 4 



^http://www. astronomy.ohio-state.edu/~sniita/xinnirsp/ 
®See the fourth item above in our list of systematics sources 
^For two observations co-added, induced features (due to 
ignored columns) should have equivalent widths of order 

integration time ratio, countrate ratio and spectral width 
of a channel, respectively 
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sometimes seen only in the longest wave- 
length range of the RGS and are a result of 
an ovcrzealous hot pixel finder. The problem 
can be ameliorated significantly by altering 
the extraction region definitions or parame- 
ters of the hot pixel finder. The problem is 
normally not noticeable in the shorter wave- 
length ranges (A < 29 A). The fluctuations 
seen in all wavelength ranges of their spec- 
trum are therefore suggestive of a different 
problem. Instead, they are probably due to 
the co-adding approach they took, which we 
referred to above. 

3. They state that pointing offsets exceeding 
15" were not included because spectral res- 
olution is degraded there. This an unneces- 
sary restriction that led to smaller counting 
statistics for them. For example, we have 
combined data from multiple offset point- 
ings, and rely on the aspect correction al- 
gorithms for the RGS to align spectral data. 
The RGS spectral resolution does not really 
degrade significantly even with ^fcw arc- 
minute scale pointing offsets. 

4. They claim that the presence of a 'detector 
feature' rendered the z ^0.011 feature un- 
measurable in the RGS. The SAS provides 
means to exclude specific detector regions, 
both from the data and from the response. 
We have shown that with proper care, sys- 
tematics arising from detector regions may 
be almost completely removed. The machin- 
ery that provides this capability is available 
and is central to the RGS branch of the SAS. 

We summarize our disagreement with Williams 
et al. (2006) as follows: While the systematics 
introduced by their method hindered spectral in- 
terpretation, those problems are not insurmount- 
able. Signatures of systematics can be distin- 
guished from instrumental PSFs, and the equiv- 
alent widths of the features reported by Nicastro 
et al. (2005b) are large enough to see by overplot- 
ting the data (c/. Fig. 4). As we have demon- 
strated above, the i? (75 data can indeed be used to 
confirm or deny presence of the weak spectral fea- 
tures that constitute the Ghandra LETGS (Nicas- 
tro et al. 2005b) detections. 



6. Conclusions 

We have analyzed a very deep spectrum of 
Mrk 421 obtained with the RGS on XMM-Newton. 
In I Ms of exposure, we collected over 26000 
counts per 50 niA resolution clement, which gives 
us high sensitivity to the detection of narrow ab- 
sorption lines. We describe the detailed procedure 
by which we reduced the data, which results in a 
continuum spectrum with noise properties that are 
dominated by statistical fluctuations. The deep 
continuum spectrum of Mrk 421 is well enough 
understood that it allows us to detect real absorp- 
tion lines of equivalent width > 2.4 mA with 99% 
confidence. 

Focusing attention on the 21.3 — 22.5 A band, 
which contains the wavelength range in which in- 
tervening O VII w absorption lines along the line 
of sight to Mrk 421 should be contained, we find 
no evidence for redshifted absorption lines. We 
place a 2 (7 upper limit of 1.9 mA on the equiv- 
alent with of any narrow absorption line over the 
quoted wavelength band. 

This finding is in clear conflict with the claim of 
the detection of narrow absorption lines at 21.85 
and 22.20 A in the deep Ghandra LETGS spec- 
trum of the source, which would naturally corre- 
spond to O VII w lines redshifted to z = 0.011 and 
z = 0.027, respectively. We find that absorption 
lines at these positions, at the equivalent widths 
claimed on the basis of the Chandra LETGS spec- 
trum, are excluded with high confidence (3.5 and 
2.8 a, respectively). The exclusion becomes only 
stronger when we use the mutually detected O VII 
z — feature as an absorption fiducial to align the 
W\ values measured across instruments as repre- 
sented in Figure 4 (5.3 and 3.8 a, respectively). 

We conclude that the detection of absorption by 
intergalactic He-like oxygen in filaments of column 
densities A'^^ « 8.5 x 10^^ cm^^ is all but excluded 
(with 99% confidence) by our data. Possible ex- 
planations for the apparent discrepancy between 
the XMM-Newton and Ghandra spectroscopic re- 
sults are investigated in an accompanying paper 
(Kaastra et al. 2006). 

The non-detection of intervening absorption, 
even in the deepest exposures that can currently 
be contemplated, on the brightest suitable extra- 
galactic continuum source, is probably not surpris- 
ing. The redshift of Mrk 421 is relatively small. 
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and the a priori probability, as predicted by recent 
cosmological gas dynamics simulations, of hav- 
ing an intervening filament with a column density 
that is readily and convincingly detectable with 
the current instrumentation, is relatively small: 
for Z - 0.1 Zq, ^ X ZMrk42i ~ 0.05 - 0.2 for 
Log{Ni) > 15.0 and 14.6, respectively (Fang et al. 
2002a; Chen et al. 2003). The most uncertain pa- 
rameter in these calculations is the absolute oxy- 
gen abundance, which can currently not be cal- 
culated with confidence from first principles. The 
null detections therefore do not yet constrain these 
calculations in a meaningful way. The detection of 
the 'missing baryons' remains a challenge, best ad- 
dressed with higher resolution, higher sensitivity 
spectrometers. 

We acknowledge several discussions with Fab- 
rizio Nicastro in which we shared our data and 
analysis, following his WHIM discovery claims. 
This work was supported by NASA. SRON is sup- 
ported financially by NWO, the Netherlands Or- 
ganization for Scientific Research. 
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Fig. 4. — A comparison of the 955 ks absorption 
spectrum toward Mrk 421 to the absorption line 
pattern of the Chandra LETGS spectrum (Nicas- 
tro et al. 2005b; Wilhams et al. 2005). 
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A. Effects of removing detector regions 

A demonstration of the iterative analysis approach is given here. The data in Figure 5 is from the dataset 
0720, and is plotted in terms of its ratio to a folded model. The top histogram is the data processed with 
automatic detection and rejection of flickering pixels, divided through with a continuum model folded through 
its corresponding response matrix. An instrumental, sharp feature located close to 21. 8 A, is identified with 
equivalent width roughly 2.5 niA. Inspection of the same data binned by hardware address reveals a single 
column with non-standard response. Because the RGS normally operates in 3x3 on-chip binning (OCB) 
mode, the charge healthy and unhealthy columns alike are added prior to digitization. In this case, it 
appears that the charge from a single unhealthy column is added to the charge of two healthy ones prior to 
signal readout. The middle panel shows the same data after reprocessing, where the channel in question was 
rejected. As a test of this method, we have also manually rejected an additional randomly selected CCD 
column whose absense is seen at a channel wavelength of ~21.37A. These data are divided by the same folded 
model and response matrix as in the top panel. The rejected column affects more than one spectral channel 
because of the chosen bin size. The bottom histogram finally shows the same data as in the middle, but 
here is divided through by the folded model using the updated response matrix, which include effects of the 
rejected columns. While the column rejection process appears to introduce some minor artifacts, substantial 
reduction in incidence of non-statistical outliers is seen. 



B. Weak spectral features arising from systematics 

We have shown above that with careful data analysis techniques, including identification of problematic 
detector locations, systematics can be reduced to a level which is comparable to the statistical limit. As 
statistics build up in very deep, large signal-to-noise observations, beating down systematics is critical in 
order for a proper identification of weak features. 

The section of the spectrum discussed in this paper is recorded on a single CCD in RGSl which, due to a 
CCD failure early in the mission, has no counterpart in RGS2. The quality of this CCD will be comparable 
to other CCD's on the RGS so the comparison between RGSl and RGS2 can give us a model independent 
estimate of the magnitude of systematic effects which may be present in RGS spectra in general and in the 
section of the spectrum analyzed in this paper in particular. 

To evaluate remaining systematic errors in RGS we have performed a model independent analysis of 
the source spectrum seen by the twin RGS instruments. By comparing the count ratio seen by the two 
instruments on a channel-by-channel basis in wavelength coordinates and comparing this to an expected 
statistical distribution, the non-statistical (or systematic) contribution may be resolved from the statistical 
population as a wing. Figure 6 displays an example of this distribution, with a Gaussian distribution (based 
on the channel counts only) overlaid to guide the eye and to help visualize the non-statistical contribution. 

Table 4 provides a tabulation of these quantities for bins of 40 mA (the approximate width of the RGS 
PSF) per CCD pair. Based on these results if one only assumes purely statistical errors, there is about a 
15% chance of falsely identifying a feature. For > 2a features this corresponds to 3 mA absorption features. 
It is stressed that the type of systematics discussed here can only result in absorption features and not in 
emission features, since all these effects can only lower the effective area and not increase it. 

Among the main causes for the non-statistical distribution in Figure 6 will be the effects of small differences 
in CTE between columns in relation to the pulse-height data selection window. Columns with clearly bad 
CTE will be manually discarded, but small CTE differences may not be recognized. They will however, 
influence the number of recorded photons due to the pulse-height selection window and hence will result in 
non-statistical count rate differences. 

Recognition of individual systematic features, which arc linked to fixed detector positions, can only be 
done with high confidence if detector coordinates change over time with respect to the wavelength coordinates 
i.e. an effective dithering of the observations. In addition, systematic features too weak to be recognized 
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Fig. 5. — A demonstration of removing detector regions in the data reduction process. A column with 
non-standard response (close to 21. SA) is rejected, along with a typical, randomly selected column (close to 
2I.37A). The three spectra depict the same data at different iterations of the reduction. 
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will be smeared-out by dithering, effectively decreasing their impact. It would therefore be advisable to 
purposely dither long RGS observations when detection of weak spectral features is desired. 



C. SAS data reduction prescriptions 

Almost all XMM observers use the standard XMM software (SAS) for their analysis. We describe in this 
section how we processed the Mrk 421 observations using the SAS and dealt with the systematic effects 
present in the data. 

The SAS does handle some aspects of CCD systematics by default, but some others require manual con- 
trol. This is only relevant in the search of weak features in spectra with a rather long observations. We used 
the SAS version 6.5, but with task 'rgsenergy' version 2.0.3 (standard SAS version 6.5 uses 'rgsenergy' ver- 
sion 2.0.2). This allows CCD-pixel dependent offsets to be subtracted (option: withdiagofFset=yes). These 
offsets are retrieved from diagnostics CCD images, sampled at regular intervals during observations. In 
addition, the default filter in the SAS to recognize and delete hot pixels was modified to delete the hot pix- 
els only, and not their neighbors (option: rejflags="BAD_SHAPE ON_BADPIX ON_WINDOW_BORDER 
BELOW.ACCEPTANCE" ). 

The SAS removes detectable, high duty cycle ("hot") pixels by default. To identify additional bad pixels 
and columns not recognized by the SAS, the data are displayed in the CCD coordinate reference frame as 
images of photon energy versus the CCD dispersion axis. Since the individual Mrk 421 observations have 
slightly different pointings, real spectral features are spread out over a range of CCD hardware coordinates, 
while CCD defects are localized in these coordinates. Such displays were used to identify local CCD defects. 
Using software written for this purpose, suspicious columns and pixels are found manually and included in 
the SAS current calibration files (CCF) to be removed during the subsequent processing. As described in the 
text, many of the columns identified in this way seem to suffer from loss of charge transfer presumably due 
to local CCD radiation damage. In the Energy versus CCD-column plot. X-ray events appear to be recorded 
at much lower energies than for the neighboring columns - resulting in a different detection efficiency into 
the local pulsehcight window. Such regions can currently only be identified manually, with the aid of high 
counting statistics observations, such as those used here. After modifying the CCF to include these bad 
columns, all data are processed again through the SAS. 

The spectra generated for each individual observation were run through the task rgsf luxer to obtain 
individually fluxed spectra. Due to the source variability these spectra were not simply added. Instead the 
spectra are all normalized over a short wavelength interval and these normalized spectra are added. The 
error is calculated taking into account the normalization. 

The final spectrum obtained in this way, using the SAS, appears to be equivalent to the spectrum obtained 
in our custom analysis. This can be seen in Figure 7. 
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Fig. 6. — The difference between RGSl and RGS2, based on background subtracted counts. The (top) 
histogram (filled circles) shows count difi'erences for single spectral bins between the two iJGS' spectrometers, 
divided by the computed statistical variance of the difference. The non-Gaussian tail component is assumed 
to be the non-statistical, systematic contribution. The overplotted continuous line, for comparison purpose, 
shows a strict Gaussian distribution, based on count statistics only. The lower distribution (filled triangles) 
show the distribution of differences between the model continuum and the actual data for the analyzed 
section of RGSl, CCD4; the CCD which is used to obtain the spectrum section discussed in this paper. 
Since this spectrum contains a spectral range of only ~1.2A, compared to the overlapping coverage between 
RGSl and RGS2 (~2lA), there is a factor of about 15 difference in number of data points binned. The 
actual distributions appear indeed to be shifted with respect to each other by a constant factor of about 10, 
confirming that our estimate of systematic errors can also be applied to the section of the spectrum discussed 
here. 
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Table 4: Consistency between RGSl and RGS2 showing only those CCD's which are operational in both 
RGSs. The flux and noise columns indicate the flux and statistical noise in the Mrk 421 spectrum for the 
specified CCD. W\ gives the equivalent width of a 1 cr feature. The last column indicates what percentage 
of the data lies outside the statistical Gaussian distribution. These results indicate that there is an overall 
probability of about 15% for a 2a feature with a typical EW of 3 mA, when the data is sampled for 40 mA 
bins. 
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Fig. 7. — Comparison between the spectrum obtained via the XMM SAS analysis package and our custom 
software. Binning is slightly different and each package identifies slightly different occasional bad pixels. 
Nevertheless, the spectra obtained are practically identical. 
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